/*******************************************************************************

				VARIANCE

  Title:        Calculate stability among files of electron/pixel.    
  Author:       Hanna Szoke
  Date:         13/09/00
  Function:     This program reads M (in range 2 - 50) sets of binary files 
		representing electron/pixels.  It calculates & writes 4 files;

		* the average at each point		(named average.bin)
		* the standard error or sqrt(variance)	(named sterror.bin)
		* error-weighted map          		(named erwm.bin)
		* an "unreliability" map		(named unrely.bin)
		
                Input parameters are read locally. UTIL.C has worker functions.

*******************************************************************************/

#include "util.h"
#include "dims.h"

static	real	*xnump ;
static	real	*average ;
static	real	*sterror ;
static	real	*erwm ;
static	real	*unrely ;

static	char	out_fn[MAXSTRING] ;

static	void	add_to_average(real *) ;  /* accumulate sum (for average)  */
static	void	add_to_sterror(real *) ;  /* accumulate sumsq (for sterror) */
static	void	allocate_var_files() ;	/* ... pointers for Np-space files */
static	void	finish_average(int) ;	/* 1/M, write average	*/
static	void	finish_sterror(int) ;	/* compute, write sterror	*/
static	void	finish_erwm() ;
static	void	finish_unrely(real) ;

void variance_main(int argc, char *argv[])
{
     void	readEpixFile(char *, real *) ;	/* Read input files */

     char	filename[MAXNFILES][MAXSTRING] ;
     int	M ;			/* number of Np file sets to handle */
     int	m ;			/* counters of files */
     real	pert_frac = 1;

  /***********************************************************
   Check for right ballpark, say hello, identify file names.
  ***********************************************************/

     M = argc - optind - 3 ;

     if (M > MAXNFILES) 
	  Fatal("Too many files!") ;

     else if (M >= 2) {
          pert_frac = atof(argv[optind+2]) ;
          for (m = 0; m < M; m++) 
               sprintf(filename[m], "%s", argv[optind+m+3]) ;
     }
     else {
	  prompt("\nPlease enter pert_frac and 2 - 8 file names: ") ;
          sscanf(terminp, "%f %s %s %s %s %s %s %s %s", &t1,
		 filename[0], filename[1], filename[2], filename[3], 
		 filename[4], filename[5], filename[6], filename[7]) ;

	  pert_frac = t1;
          M = m = 0 ;
          while ((int)strlen(filename[m]) > 0)   {
	       M++; 
	       m++;
          }
          if (M < 2) 
	       Fatal("Too few files!") ;

          sprintf(message, "Running %s with pert_frac %f on %d files", 
	          caller, pert_frac, M) ; 
          printTwice(message) ;
     }

     hello(caller) ;

	/************************************************************
	Fetch input conditions, set up arrays in physical space.
	Read each binary el/voxel file, add to an average array.
	************************************************************/

     readBasicInput() ;

     allocate_var_files() ;

     if (HighRes) {
	  printTwice("\n\t\tNOTE:") ;
	  printTwice(
     "Variance does not currently handle lists of high-resolution points.\n") ;
     }

     sprintf(message, "\nNumber of files = %d, Pert_frac = %g\n", 
             M, pert_frac);
     printTwice(message) ;

     printTwice("Calculating average among input files.\n") ;

     for (m = 0; m < M; m++)  {
          readEpixFile(filename[m], xnump) ;
          add_to_average(xnump) ;
     }
     finish_average(M) ;

	/************************************************************
	Compute and write standard error; this requires that we 
	reread the original files.
	************************************************************/

     printTwice("Calculating standard error among input files.\n") ;

     for (m = 0; m < M; m++)  {
          readEpixFile(filename[m], xnump) ;
          add_to_sterror(xnump) ;
     }
     finish_sterror(M) ;

	/************************************************************
	Compute and write the error-weighted map and "unreliability".
	************************************************************/

     printTwice(
     "Calculating the error-weighted average and unreliabiility maps.\n") ;

     finish_erwm() ;
     finish_unrely(pert_frac) ;

     return ;
}

void	allocate_var_files()
{

     xnump = (real *) e_malloc(Nptotal*sizeof(real), caller) ;
     average = (real *) e_malloc(Nptotal*sizeof(real), caller) ;
     sterror = (real *) e_malloc(Nptotal*sizeof(real), caller) ;
     erwm = (real *) e_malloc(Nptotal*sizeof(real), caller) ;
     unrely = (real *) e_malloc(Nptotal*sizeof(real), caller) ;

}
void	add_to_average(real *nump)		/* accumulate nump */  
{
     int 	n ;
     static	int	first = TRUE ;

     real	*pnump, *pave ;

     if (first) {
          first = FALSE ;
          for (n = 0, pave = average; n <Nptotal; n++, pave++)
               *pave = 0.;
     }

     for (n = 0, pave = average, pnump = nump; n <Nptotal; n++, pave++,pnump++)
          *pave += *pnump ;

     return;
}

void	finish_average(int M)
{
     void       writeEpixFile(char *, real *) ;

     int 	n ;
     real	*pave ;
 
     for (n = 0, pave = average; n <Nptotal; n++, pave++)
          *pave /= M ;

     sprintf(out_fn, "average") ;
     writeEpixFile(out_fn, average) ;

     return;
}

void	add_to_sterror(real *nump)		/* accumulate nump * nump */  
{
     int 	n ;
     static	int	first = TRUE ;

     real	*pnump, *pste ;

     if (first) {
          first = FALSE ;
          for (n = 0, pste = sterror; n <Nptotal; n++, pste++)
               *pste = 0. ;
     }

     for (n = 0, pste = sterror, pnump = nump; n <Nptotal; n++, pste++,pnump++)
          *pste += *pnump * *pnump ;

     return;
}

void	finish_sterror(int M)
{
     void       writeEpixFile() ;

     int 	n ;
     real	*pave, *pste ;
 
     for (n = 0, pste = sterror, pave = average; n <Nptotal; 
          n++, pste++, pave++) {

          *pste -= M * *pave * *pave ;
          *pste /= (M-1) ;
          *pste = sqrt(fabs(*pste)) ;
     }

     sprintf(out_fn, "sterror") ;
     writeEpixFile(out_fn, sterror) ;

     return;
}

void	finish_erwm()
{
     void       writeEpixFile() ;

     int	n ;
     real	*pave, *pste, *per ;

     for (n = 0, pste = sterror, pave = average, per = erwm; n <Nptotal; 
          n++, pste++, pave++, per++) {

          if (*pste > 1.e-5)
              *per = (*pave * *pave) / (*pave + *pste) ;
          else
              *per = *pave ;
     }
     sprintf(out_fn, "erwm") ;
     writeEpixFile(out_fn, erwm) ;

     return;
}

void	finish_unrely(real pert_frac)
{
     void       writeEpixFile() ;

     int	n ;
     real	*pave, *pste, *pun ;

     for (n = 0, pste = sterror, pave = average, pun = unrely; n <Nptotal; 
          n++, pste++, pave++, pun++) 

         if (*pste > 1.e-5)
              *pun = *pste / (pert_frac * *pave + *pste) ;
	 else
              *pun = 0 ;
						
     sprintf(out_fn, "unrely") ;
     writeEpixFile(out_fn, unrely) ;

     return;
}
